Code to reproduce results of the manuscript ‘Kidney Failure Prediction: Multicenter External Validation of the KFRE Model in Patients with CKD Stages 3-4 in Peru’
Author
Percy Soto-Becerra
Published
March 29, 2023
1 Introduction
This document presents the code and results of the main analysis shown in the article.
2 Setup
Mostrar código
rm(list =ls())# Use pacman to check whether packages are installed, if not loadif(!require("pacman"))install.packages("pacman")library(pacman)# Unload all package to begin in a session with only base packagespacman::p_unload("all")# Install packagespacman::p_load(here, skimr, survival,rms,cmprsk,riskRegression,mstate,pseudo,pec,riskRegression,plotrix,knitr,splines,kableExtra,flextable,gtsummary,boot,tidyverse,rsample,gridExtra,webshot, patchwork,survival, cmprsk, survminer, ggsci, cowplot, scales, patchwork, labelled, glue, dcurves, broom, downlit, xml2, gghalves, devtools, htmltools, gghalves, ggtext, DiagrammeR, gt, janitor)if(!require("smplot2"))devtools::install_github('smin95/smplot2', force =TRUE)# Import datadata<-readRDS(here::here("Data/Derived/data_derived.rds"))# Subset patients with CKD Stages 3a-3b-4data%>%filter(ckd_stage=="Stages 3-4")->dataA# Subset patients with CKD Stages 3b-4data%>%filter(ckd_stage2=="Stages 3b-4")->dataB
3 Descriptive analysis
In total, 7519 patients were included in the analysis of the group with stages 3a-4 CKD, and 2798 in the group with stages 3b-4. All patients had outcome data. The number of events was greater than 100 in all populations and outcomes, except for renal failure at 2 years in the population with stages 3b-4 (n = 88), so estimates in this group should be interpreted with more caution. Specifically, in the 3a-3b-4 group, 114 (1.5%) patients developed renal failure at 2 years and 239 (3.2%) at 5 years. Many more patients died without experiencing renal failure: 563 (7.5%) at 2 years and 1400 (18.6%) at 5 years.
Regarding the group restricted to stages 3b-4, 88 (3.1%) patients developed renal failure at 2 years and 182 (6.5%) at 5 years. Similarly, many more patients died without experiencing renal failure: 300 (10.7%) at 2 years and 683 (24.4%) at 5 years.
The median observation time was 4.9 years and the maximum follow-up was 7.8 years in the 3a-4 group. The median and maximum observation times were similar in the 3b-4 group, at 4.9 and 7.8 years, respectively.
3.1 Fig 1
Mostrar código
# Create grid of 100 x 100----data_flow<-tibble(x =1:100, y =1:100)data_flow%>%ggplot(aes(x, y))+scale_x_continuous(minor_breaks =seq(1, 100, 1))+scale_y_continuous(minor_breaks =seq(1, 100, 1))+theme_linedraw()->p# Create boxes----# box_xmin<-33-20box_xmax<-75-20box_ymin<-94box_ymax<-99box_size<-0.25text_param<-function(box_min, box_max){mean(c(box_min, box_max))}text_size<-2.5p+# Col 0----## Level 1----geom_rect(xmin =box_xmin, xmax =box_xmax, ymin =box_ymin-1, ymax =box_ymax+1, color ="black", fill ="white", size =box_size)+annotate('text', x =text_param(box_xmin, box_xmax), y =text_param(box_ymin-1, box_ymax+1), label ='Total patients in VISARE database\n(n = 57,308)', size =text_size)+## Level 2----geom_rect(xmin =box_xmin+5, xmax =box_xmax-5, ymin =box_ymin-41, ymax =box_ymax-39, color ="black", fill ="white", size =box_size)+annotate('text', x =text_param(box_xmin+5, box_xmax-5), y =text_param(box_ymin-41, box_ymax-39), label ='Total of patients with CKD\n(n = 22,744)', size =text_size)+# Col -1----geom_rect(xmin =box_xmin-13, xmax =box_xmax-27, ymin =box_ymin-81, ymax =box_ymax-79, color ="black", fill ="white", size =box_size)+annotate('text', x =text_param(box_xmin-13, box_xmax-27), y =text_param(box_ymin-81, box_ymax-79), label ='Patients with CKD 3a-4\n(n = 7,519)', size =text_size)+# Col +1----## Level 1----geom_rect(xmin =box_xmin+23, xmax =box_xmax+47, ymin =box_ymin-19, ymax =box_ymax-12, color ="black", fill ="white", size =box_size)+annotate('text', x =text_param(box_xmin+23, box_xmax+47), y =text_param(box_ymin-19, box_ymax-12), label ='34,564 screened patients did not CKD in any stage', size =text_size)+## Level 2----geom_rect(xmin =box_xmin+23, xmax =box_xmax+47, ymin =box_ymin-19-45, ymax =box_ymax-12-45, color ="black", fill ="white", size =box_size)+annotate('text', x =text_param(box_xmin+23, box_xmax+47), y =text_param(box_ymin-19-45, box_ymax-12-45), label ='Not elegible (n = 15,225) \n 8,426 lacks of data on urine albumin and urine creatinin \n 6,799 were in other CKD stages (1, 2, or 5)', size =text_size)+## Level 3----geom_rect(xmin =box_xmin+27, xmax =box_xmax+13, ymin =box_ymin-81 , ymax =box_ymax-79, color ="black", fill ="white", size =box_size)+annotate('text', x =text_param(box_xmin+27, box_xmax+13), y =text_param(box_ymin-81, box_ymax-79), label ='Patients with CKD 3b-4\n(n = 2,798)', size =text_size)+# vertical arrowgeom_segment(x =text_param(box_xmin, box_xmax), xend =text_param(box_xmin, box_xmax), y =box_ymin-1, yend =box_ymax-39, size =0.15, linejoin ="mitre", lineend ="butt", arrow =arrow(length =unit(1, "mm"), type ="closed"))+geom_segment(x =text_param(box_xmin, box_xmax), xend =text_param(box_xmin, box_xmax), y =box_ymin-41, yend =25, size =0.15, linejoin ="mitre", lineend ="butt")+# horizontal arrow 1-->geom_segment(x =text_param(box_xmin, box_xmax), xend =box_xmin+23, y =text_param(box_ymin-19, box_ymax-12), yend =text_param(box_ymin-19, box_ymax-12), size =0.15, linejoin ="mitre", lineend ="butt", arrow =arrow(length =unit(1, "mm"), type ="closed"))+# horizontal arrow 2-->geom_segment(x =text_param(box_xmin, box_xmax), xend =box_xmin+23, y =text_param(box_ymin-19-45, box_ymax-12-45), yend =text_param(box_ymin-19-45, box_ymax-12-45), size =0.15, linejoin ="mitre", lineend ="butt", arrow =arrow(length =unit(1, "mm"), type ="closed"))+# horizontal segment --geom_segment(x =text_param(box_xmin-13, box_xmax-27), xend =text_param(box_xmin+27, box_xmax+13), y =25, yend =25, size =0.15, linejoin ="mitre", lineend ="butt")+# vertical arrow -->geom_segment(x =text_param(box_xmin-13, box_xmax-27), xend =text_param(box_xmin-13, box_xmax-27), y =25, yend =box_ymax-79, size =0.15, linejoin ="mitre", lineend ="butt", arrow =arrow(length =unit(1, "mm"), type ="closed"))+# vertical arrow -->geom_segment(x =text_param(box_xmin+27, box_xmax+13), xend =text_param(box_xmin+27, box_xmax+13), y =25, yend =box_ymax-79, size =0.15, linejoin ="mitre", lineend ="butt", arrow =arrow(length =unit(1, "mm"), type ="closed"))+theme_void()+theme(plot.background =element_rect(fill ="white"))->plot_flowchartggsave(filename ="plot_flowchart.png", plot =plot_flowchart, device ="png", path =here("Figures"), scale =1, width =12, height =12, units ="cm", dpi =600)
Table 1 provides a description of the study population. Table S1 (see section “tableS1”) details the characteristics of the study population with stages 3a-3b-4 according to the occurrence of the outcome of interest, renal failure at 2 and 5 years. The number of events of interest at 2 years was reported, and Table S2 (see section “tableS2”) shows the same information for the 3b-4 population. Table S3 (see section “tableS3”) describes the characteristics according to stages 3a, 3b, and 4 separately. With this stratification, the numbers of cases of renal failure at 2 years were very low for the subpopulations with stages 3a (n = 26), 3b (n = 36), and 4 (n = 52). Likewise, the numbers of cases at 5 years were very low for the subpopulations with stages 3a (n = 57), 3b (n = 81), and 4 (n = 101). Therefore, it was not reliable to perform predictive performance evaluation in these subgroups.
Table 1. Baseline characteristics of the study population according CKD Stages
CKD Stages 3a-3b-4
CKD Stages 3b-4
Characteristic
N = 7,519
N = 2,798
Sex
Female
4,107 (54.6%)
1,398 (50.0%)
Male
3,412 (45.4%)
1,400 (50.0%)
Age (years)
Mean (SD)
74.0 (10.2)
75.6 (10.6)
Median (IQR)
75.0 (68.0 - 81.0)
77.0 (70.0 - 83.0)
Range
23.0 - 97.0
23.0 - 97.0
Hypertension
4,486 (59.7%)
1,636 (58.5%)
Diabetes Mellitus
1,845 (24.5%)
674 (24.1%)
Persistent albuminuria categories
A1
4,772 (63.5%)
1,494 (53.4%)
A2
2,018 (26.8%)
860 (30.7%)
A3
729 (9.7%)
444 (15.9%)
GFR categories
G3a
4,721 (62.8%)
G3b
2,207 (29.4%)
2,207 (78.9%)
G4
591 (7.9%)
591 (21.1%)
CKD KDIGO classification
Moderately increased risk
3,278 (43.6%)
High risk
2,460 (32.7%)
1,302 (46.5%)
Very high risk
1,781 (23.7%)
1,496 (53.5%)
Serum Creatinine (mg/dL)
Mean (SD)
1.4 (0.4)
1.7 (0.5)
Median (IQR)
1.3 (1.1 - 1.5)
1.6 (1.4 - 1.9)
Range
0.8 - 3.9
1.1 - 3.9
eGFR using CKD-EPI, ml/min/1.73m2
Mean (SD)
46.2 (9.8)
35.7 (7.3)
Median (IQR)
48.7 (40.4 - 53.8)
37.3 (31.4 - 41.7)
Range
15.0 - 60.0
15.0 - 45.0
Albumin-to-creatinine ratio, mg/g
Mean (SD)
248.6 (3,044.4)
334.1 (3,050.8)
Median (IQR)
14.6 (4.5 - 66.1)
26.0 (6.5 - 153.8)
Range
0.0 - 144,870.6
0.0 - 144,870.6
Urine Albumin (mg/ml)
Mean (SD)
8.3 (28.3)
13.1 (36.3)
Median (IQR)
1.0 (0.3 - 4.0)
1.7 (0.4 - 10.1)
Range
0.0 - 658.0
0.0 - 658.0
Urine Creatinine (mg/dl)
Mean (SD)
72.4 (47.5)
71.4 (47.0)
Median (IQR)
63.3 (41.3 - 86.0)
64.9 (43.3 - 85.0)
Range
0.1 - 722.1
0.7 - 722.1
Death at 2 years
640 (8.5%)
358 (12.8%)
Outcome at 2 years
Alive w/o Kidney Failure
6,842 (91.0%)
2,410 (86.1%)
Death w/o Kidney Failure
563 (7.5%)
300 (10.7%)
Kidney Failure
114 (1.5%)
88 (3.1%)
Death at 5 years
1,539 (20.5%)
784 (28.0%)
Outcome at 5 years
Alive w/o Kidney Failure
5,880 (78.2%)
1,933 (69.1%)
Death w/o Kidney Failure
1,400 (18.6%)
683 (24.4%)
Kidney Failure
239 (3.2%)
182 (6.5%)
1 n (%)
4 Cumulative incidence function for competing risks data
Figure Figure 1 shows the cumulative incidence curves of renal failure and pre-renal failure death in the study patients.
Mostrar código
# Seleccion del grupo 3a-4----cuminc(ftime =dataA$time5y, fstatus =dataA$eventd5ylab, cencode ="Alive w/o Kidney Failure")->cifciplotdat<-cif%>%list_modify("Tests"=NULL)%>%map_df(`[`, c("time", "est", "var"), .id ="id")%>%mutate(id =recode(id, "1 Death w/o Kidney Failure"="Death w/o Kidney Failure", "1 Kidney Failure"="Kidney Failure"), ll =est-1.96*sqrt(var), ul =est+1.96*sqrt(var))%>%rename( event =id)ciplotdat%>%ggplot(aes(x =time, y =est))+geom_ribbon(aes(ymin =ll, ymax =ul, fill =event), alpha =0.25, linetype =0)+geom_step(lwd =0.5, aes(color =event))+theme_survminer()+scale_y_continuous(labels =scales::percent, limits =c(0, 1))+labs(x ="Years", y ="Cumulative incidence", title ="A CKD Stages 3a-3b-4")+theme(legend.position ="top", legend.title =element_blank(), legend.background =element_rect(fill ="white"), legend.key.size =unit(0.2, "cm"))->g1ciplotdat%>%filter(event=="Kidney Failure")%>%ggplot(aes(x =time, y =est))+geom_ribbon(aes(ymin =ll, ymax =ul), fill ="#89E1E3", alpha =0.1, linetype =0)+geom_step(lwd =0.5, color ="#89E1E3")+theme_survminer()+ylim(c(0, 0.025))+scale_y_continuous(labels =scales::percent, limits =c(0, 0.04))+labs(x ="", y ="", title ="")->g2kf_fit<-survfit(Surv(time5y, ifelse(eventd5y!=0, 1, 0))~1, data =dataA)num<-ggsurvplot( fit =kf_fit, risk.table ="nrisk_cumevents", risk.table.y.text =FALSE, risk.table.y.text.col =FALSE, tables.y.text =FALSE, tables.y.text.col =FALSE, ylab ="Years", risk.table.fontsize =3.2, tables.theme =theme_survminer(font.main =10))cowplot::plot_grid(g1, num$table+theme_cleantable(), nrow =2, rel_heights =c(4, 1), align ="v", axis ="b")->g3g3+inset_element(g2, 0.15, 0.43, 1, 0.856, align_to ='full', ignore_tag =TRUE)->plot_cif_mh# Seleccion del grupo 3b-4----cuminc(ftime =dataB$time5y, fstatus =dataB$eventd5ylab, cencode ="Alive w/o Kidney Failure")->cifciplotdat<-cif%>%list_modify("Tests"=NULL)%>%map_df(`[`, c("time", "est", "var"), .id ="id")%>%mutate(id =recode(id, "1 Death w/o Kidney Failure"="Death w/o Kidney Failure", "1 Kidney Failure"="Kidney Failure"), ll =est-1.96*sqrt(var), ul =est+1.96*sqrt(var))%>%rename( event =id)ciplotdat%>%ggplot(aes(x =time, y =est))+geom_ribbon(aes(ymin =ll, ymax =ul, fill =event), alpha =0.25, linetype =0)+geom_step(lwd =0.5, aes(color =event))+theme_survminer()+scale_y_continuous(labels =scales::percent, limits =c(0, 1))+labs(x ="Years", y ="Cumulative incidence", title ="B CKD Stages 3b-4")+theme(legend.position ="top", legend.title =element_blank(), legend.background =element_rect(fill ="white"), legend.key.size =unit(0.2, "cm"))->g1ciplotdat%>%filter(event=="Kidney Failure")%>%ggplot(aes(x =time, y =est))+geom_ribbon(aes(ymin =ll, ymax =ul), fill ="#89E1E3", alpha =0.1, linetype =0)+geom_step(lwd =0.5, color ="#89E1E3")+theme_survminer()+ylim(c(0, 0.15))+scale_y_continuous(labels =scales::percent, limits =c(0, 0.08))+labs(x ="", y ="", title ="")->g2kf_fit<-survfit(Surv(time5y, ifelse(eventd5y!=0, 1, 0))~1, data =dataB)num<-ggsurvplot( fit =kf_fit, risk.table ="nrisk_cumevents", risk.table.y.text =FALSE, risk.table.y.text.col =FALSE, tables.y.text =FALSE, tables.y.text.col =FALSE, ylab ="Years", risk.table.fontsize =3.2, tables.theme =theme_survminer(font.main =10))cowplot::plot_grid(g1, num$table+theme_cleantable(), nrow =2, rel_heights =c(4, 1), align ="v", axis ="b")->g3g3+inset_element(g2, 0.15, 0.51, 1, 0.856, align_to ='full', ignore_tag =TRUE)->plot_cif_vh(plot_cif_mh/plot_cif_vh)+plot_annotation(tag_levels ='A')->plot_cifggsave(filename ="Plot_CIF.png", plot =plot_cif, device ="png", path =here("Figures"), dpi =300, scale =2, width =8.5, height =14, units ="cm", bg ="white")
4.1 Fig 2
Mostrar código
plot_cif
Figure 1: Cumulative incidence curves of the observed outcome probabilities in the population study for kidney failure or deathd without kidney failure
5 Predictive Performance
5.1 Calibration
Mean calibration: OE ratio
Mostrar código
# Calibration (O/E) -------------------------------------------------------# Seleccion del grupo: Stages 3-4----vdata<-dataA%>%select(id, risk2y, risk5y, time5y, eventd5y, time, eventd)%>%drop_na()%>%mutate(status =factor(eventd5y, levels =c(0, 1, 2), labels =c("Alive w/o Kidney Failure", "Kidney Failure", "Death w/o Kidney Failure")))primary_event<-1# A 2 años----horizon<-2vdata$pred<-vdata$risk2y# First calculate Aalen-Johansen estimate (as 'observed')obj<-summary(survfit(Surv(time5y, status)~1, data =vdata), times =horizon)aj<-list("obs"=obj$pstate[, primary_event+1],"se"=obj$std.err[, primary_event+1])# Calculate O/EOE<-aj$obs/mean(vdata$pred)# For the confidence interval we use method proposed in Debray et al. (2017) doi:10.1136/bmj.i6460k<-2alpha<-0.05OE_summary<-cbind("OE"=OE,"Lower .95"=exp(log(OE-qnorm(1-alpha/2)*aj$se/aj$obs)),"Upper .95"=exp(log(OE+qnorm(1-alpha/2)*aj$se/aj$obs)))OE_summary2a<-round(OE_summary, k)avg_obs<-cbind("avg_obs"=aj$obs*100,"Lower .95"=100*(aj$obs-qnorm(1-alpha/2)*aj$se),"Upper .95"=100*(aj$obs+qnorm(1-alpha/2)*aj$se))avg_pred<-cbind("avg_pred"=mean(vdata$pred)*100)avg_obs2a<-round(avg_obs, k)avg_pred2a<-round(avg_pred, k)# A 5 años----horizon<-5# Add estimated risk and complementary log-log of it to datasetvdata$pred<-vdata$risk5y# First calculate Aalen-Johansen estimate (as 'observed')obj<-summary(survfit(Surv(time5y, status)~1, data =vdata), times =horizon)aj<-list("obs"=obj$pstate[, primary_event+1],"se"=obj$std.err[, primary_event+1])# Calculate O/EOE<-aj$obs/mean(vdata$pred)# For the confidence interval we use method proposed in Debray et al. (2017) doi:10.1136/bmj.i6460k<-2alpha<-0.05OE_summary<-cbind("OE"=OE,"Lower .95"=exp(log(OE-qnorm(1-alpha/2)*aj$se/aj$obs)),"Upper .95"=exp(log(OE+qnorm(1-alpha/2)*aj$se/aj$obs)))OE_summary5a<-round(OE_summary, k)avg_obs<-cbind("avg_obs"=aj$obs*100,"Lower .95"=100*(aj$obs-qnorm(1-alpha/2)*aj$se),"Upper .95"=100*(aj$obs+qnorm(1-alpha/2)*aj$se))avg_pred<-cbind("avg_pred"=mean(vdata$pred)*100)avg_obs5a<-round(avg_obs, k)avg_pred5a<-round(avg_pred, k)# Seleccion del grupo: Stages 3b-4----vdata<-dataB%>%select(id, risk2y, risk5y, time5y, eventd5y, eventd, time)%>%drop_na()%>%mutate(status =factor(eventd5y, levels =c(0, 1, 2), labels =c("Alive w/o Kidney Failure", "Kidney Failure", "Death w/o Kidney Failure")))primary_event<-1# A 2 años----horizon<-2# Add estimated risk and complementary log-log of it to datasetvdata$pred<-vdata$risk2y# First calculate Aalen-Johansen estimate (as 'observed')obj<-summary(survfit(Surv(time5y, status)~1, data =vdata), times =horizon)aj<-list("obs"=obj$pstate[, primary_event+1],"se"=obj$std.err[, primary_event+1])# Calculate O/EOE<-aj$obs/mean(vdata$pred)# For the confidence interval we use method proposed in Debray et al. (2017) doi:10.1136/bmj.i6460k<-2alpha<-0.05OE_summary<-cbind("OE"=OE,"Lower .95"=exp(log(OE-qnorm(1-alpha/2)*aj$se/aj$obs)),"Upper .95"=exp(log(OE+qnorm(1-alpha/2)*aj$se/aj$obs)))OE_summary2b<-round(OE_summary, k)avg_obs<-cbind("avg_obs"=aj$obs*100,"Lower .95"=100*(aj$obs-qnorm(1-alpha/2)*aj$se),"Upper .95"=100*(aj$obs+qnorm(1-alpha/2)*aj$se))avg_pred<-cbind("avg_pred"=mean(vdata$pred)*100)avg_obs2b<-round(avg_obs, k)avg_pred2b<-round(avg_pred, k)# A 5 años----horizon<-5# Add estimated risk and complementary log-log of it to datasetvdata$pred<-vdata$risk5y# First calculate Aalen-Johansen estimate (as 'observed')obj<-summary(survfit(Surv(time5y, status)~1, data =vdata), times =horizon)aj<-list("obs"=obj$pstate[, primary_event+1],"se"=obj$std.err[, primary_event+1])# Calculate O/EOE<-aj$obs/mean(vdata$pred)# For the confidence interval we use method proposed in Debray et al. (2017) doi:10.1136/bmj.i6460k<-2alpha<-0.05OE_summary<-cbind("OE"=OE,"Lower .95"=exp(log(OE-qnorm(1-alpha/2)*aj$se/aj$obs)),"Upper .95"=exp(log(OE+qnorm(1-alpha/2)*aj$se/aj$obs)))OE_summary5b<-round(OE_summary, k)avg_obs<-cbind("avg_obs"=aj$obs*100,"Lower .95"=100*(aj$obs-qnorm(1-alpha/2)*aj$se),"Upper .95"=100*(aj$obs+qnorm(1-alpha/2)*aj$se))avg_pred<-cbind("avg_pred"=mean(vdata$pred)*100)avg_obs5b<-round(avg_obs, k)avg_pred5b<-round(avg_pred, k)
Weak calibration: Calibration intercept and Calibration slope
Mostrar código
# Seleccion del grupo: Stages 3-4----vdata<-dataA%>%select(id, risk2y, risk5y, time5y, eventd5y, eventd, time)%>%drop_na()%>%mutate(status =factor(eventd5y, levels =c(0, 1, 2), labels =c("Alive w/o Kidney Failure", "Kidney Failure", "Death w/o Kidney Failure")))primary_event<-1# A 2 años----horizon<-2# Add estimated risk and complementary log-log of it to datasetvdata$pred<-vdata$risk2y# First compute riskRegression::Score()score_vdata<-Score(list("csh_validation"=vdata$pred), formula =Hist(time, eventd)~1, cens.model ="km", data =vdata, conf.int =TRUE, times =horizon,# metrics = c("auc", "brier"), summary =c("ipa"), cause =primary_event, plots ="calibration")# Use pseudo-observations calculated by Score() (can alternatively use pseudo::pseudoci)pseudos<-data.frame(score_vdata$Calibration$plotframe)# Note:# - 'pseudos' is the data.frame with ACTUAL pseudo-observations, not the smoothed ones# - Column ID is not the id in vdata; it is just a number assigned to each row of# the original validation data sorted by time and event indicatorpseudos$cll_pred<-log(-log(1-pseudos$risk))# add the cloglog risk ests# Fit model for calibration interceptfit_cal_int<-geese(pseudovalue~offset(cll_pred), data =pseudos, id =ID, scale.fix =TRUE, family =gaussian, mean.link ="cloglog", corstr ="independence", jack =TRUE)# Fit model for calibration slopefit_cal_slope<-geese(pseudovalue~offset(cll_pred)+cll_pred, data =pseudos, id =ID, scale.fix =TRUE, family =gaussian, mean.link ="cloglog", corstr ="independence", jack =TRUE)# Perform joint test on intercept and slopebetas<-fit_cal_slope$betavcov_mat<-fit_cal_slope$vbetawald<-drop(betas%*%solve(vcov_mat)%*%betas)# pchisq(wald, df = 2, lower.tail = FALSE)k<-2res_cal<-rbind(# Value, confidence interval and test for calibration intercept"Intercept, Stages 3-4, 2 year"=with(summary(fit_cal_int)$mean,c("estimate"=estimate, `2.5 %` =estimate-qnorm(0.975)*san.se, `97.5 %` =estimate+qnorm(0.975)*san.se)),# Value, confidence interval and test for calibration slope"Slope, Stages 3-4, 2 year"=with(summary(fit_cal_slope)$mean["cll_pred", ],c("estimate"=1+estimate, `2.5 %` =1+(estimate-qnorm(0.975)*san.se), `97.5 %` =1+(estimate+qnorm(0.975)*san.se))))res_cal2a<-round(res_cal, k)# A 5 años----horizon<-5# Add estimated risk and complementary log-log of it to datasetvdata$pred<-vdata$risk5y# First compute riskRegression::Score()score_vdata<-Score(list("csh_validation"=vdata$pred), formula =Hist(time, eventd)~1, cens.model ="km", data =vdata, conf.int =TRUE, times =horizon,# metrics = c("auc", "brier"), summary =c("ipa"), cause =primary_event, plots ="calibration")# Use pseudo-observations calculated by Score() (can alternatively use pseudo::pseudoci)pseudos<-data.frame(score_vdata$Calibration$plotframe)# Note:# - 'pseudos' is the data.frame with ACTUAL pseudo-observations, not the smoothed ones# - Column ID is not the id in vdata; it is just a number assigned to each row of# the original validation data sorted by time and event indicatorpseudos$cll_pred<-log(-log(1-pseudos$risk))# add the cloglog risk ests# Fit model for calibration interceptfit_cal_int<-geese(pseudovalue~offset(cll_pred), data =pseudos, id =ID, scale.fix =TRUE, family =gaussian, mean.link ="cloglog", corstr ="independence", jack =TRUE)# Fit model for calibration slopefit_cal_slope<-geese(pseudovalue~offset(cll_pred)+cll_pred, data =pseudos, id =ID, scale.fix =TRUE, family =gaussian, mean.link ="cloglog", corstr ="independence", jack =TRUE)# Perform joint test on intercept and slopebetas<-fit_cal_slope$betavcov_mat<-fit_cal_slope$vbetawald<-drop(betas%*%solve(vcov_mat)%*%betas)# pchisq(wald, df = 2, lower.tail = FALSE)k<-2res_cal<-rbind(# Value, confidence interval and test for calibration intercept"Intercept, Stages 3-4, 5 year"=with(summary(fit_cal_int)$mean,c("estimate"=estimate, `2.5 %` =estimate-qnorm(0.975)*san.se, `97.5 %` =estimate+qnorm(0.975)*san.se)),# Value, confidence interval and test for calibration slope"Slope, Stages 3-4, 5 year"=with(summary(fit_cal_slope)$mean["cll_pred", ],c("estimate"=1+estimate, `2.5 %` =1+(estimate-qnorm(0.975)*san.se), `97.5 %` =1+(estimate+qnorm(0.975)*san.se))))res_cal5a<-round(res_cal, k)# Seleccion del grupo: Stages 3b-4----vdata<-dataB%>%select(id, risk2y, risk5y, time5y, eventd5y, time, eventd)%>%drop_na()%>%mutate(status =factor(eventd5y, levels =c(0, 1, 2), labels =c("Alive w/o Kidney Failure", "Kidney Failure", "Death w/o Kidney Failure")))primary_event<-1# A 2 años----horizon<-2# Add estimated risk and complementary log-log of it to datasetvdata$pred<-vdata$risk2y# First compute riskRegression::Score()score_vdata<-Score(list("csh_validation"=vdata$pred), formula =Hist(time, eventd)~1, cens.model ="km", data =vdata, conf.int =TRUE, times =horizon,# metrics = c("auc", "brier"), summary =c("ipa"), cause =primary_event, plots ="calibration")# Use pseudo-observations calculated by Score() (can alternatively use pseudo::pseudoci)pseudos<-data.frame(score_vdata$Calibration$plotframe)# Note:# - 'pseudos' is the data.frame with ACTUAL pseudo-observations, not the smoothed ones# - Column ID is not the id in vdata; it is just a number assigned to each row of# the original validation data sorted by time and event indicatorpseudos$cll_pred<-log(-log(1-pseudos$risk))# add the cloglog risk ests# Fit model for calibration interceptfit_cal_int<-geese(pseudovalue~offset(cll_pred), data =pseudos, id =ID, scale.fix =TRUE, family =gaussian, mean.link ="cloglog", corstr ="independence", jack =TRUE)# Fit model for calibration slopefit_cal_slope<-geese(pseudovalue~offset(cll_pred)+cll_pred, data =pseudos, id =ID, scale.fix =TRUE, family =gaussian, mean.link ="cloglog", corstr ="independence", jack =TRUE)# Perform joint test on intercept and slopebetas<-fit_cal_slope$betavcov_mat<-fit_cal_slope$vbetawald<-drop(betas%*%solve(vcov_mat)%*%betas)# pchisq(wald, df = 2, lower.tail = FALSE)k<-2res_cal<-rbind(# Value, confidence interval and test for calibration intercept"Intercept, Stages 3b-4, 2 year"=with(summary(fit_cal_int)$mean,c("estimate"=estimate, `2.5 %` =estimate-qnorm(0.975)*san.se, `97.5 %` =estimate+qnorm(0.975)*san.se)),# Value, confidence interval and test for calibration slope"Slope, Stages 3b-4, 2 year"=with(summary(fit_cal_slope)$mean["cll_pred", ],c("estimate"=1+estimate, `2.5 %` =1+(estimate-qnorm(0.975)*san.se), `97.5 %` =1+(estimate+qnorm(0.975)*san.se))))res_cal2b<-round(res_cal, k)# A 5 años----horizon<-5vdata$pred<-vdata$risk5y# First compute riskRegression::Score()score_vdata<-Score(list("csh_validation"=vdata$pred), formula =Hist(time, eventd)~1, cens.model ="km", data =vdata, conf.int =TRUE, times =horizon,# metrics = c("auc", "brier"), summary =c("ipa"), cause =primary_event, plots ="calibration")# Use pseudo-observations calculated by Score() (can alternatively use pseudo::pseudoci)pseudos<-data.frame(score_vdata$Calibration$plotframe)# Note:# - 'pseudos' is the data.frame with ACTUAL pseudo-observations, not the smoothed ones# - Column ID is not the id in vdata; it is just a number assigned to each row of# the original validation data sorted by time and event indicatorpseudos$cll_pred<-log(-log(1-pseudos$risk))# add the cloglog risk ests# Fit model for calibration interceptfit_cal_int<-geese(pseudovalue~offset(cll_pred), data =pseudos, id =ID, scale.fix =TRUE, family =gaussian, mean.link ="cloglog", corstr ="independence", jack =TRUE)# Fit model for calibration slopefit_cal_slope<-geese(pseudovalue~offset(cll_pred)+cll_pred, data =pseudos, id =ID, scale.fix =TRUE, family =gaussian, mean.link ="cloglog", corstr ="independence", jack =TRUE)# Perform joint test on intercept and slopebetas<-fit_cal_slope$betavcov_mat<-fit_cal_slope$vbetawald<-drop(betas%*%solve(vcov_mat)%*%betas)# pchisq(wald, df = 2, lower.tail = FALSE)k<-2res_cal<-rbind(# Value, confidence interval and test for calibration intercept"Intercept, Stages 3b-4, 5 year"=with(summary(fit_cal_int)$mean,c("estimate"=estimate, `2.5 %` =estimate-qnorm(0.975)*san.se, `97.5 %` =estimate+qnorm(0.975)*san.se)),# Value, confidence interval and test for calibration slope"Slope, Stages 3b-4, 5 year"=with(summary(fit_cal_slope)$mean["cll_pred", ],c("estimate"=1+estimate, `2.5 %` =1+(estimate-qnorm(0.975)*san.se), `97.5 %` =1+(estimate+qnorm(0.975)*san.se))))res_cal5b<-round(res_cal, k)
Calibration-in-the-large
OE ratio is a measure of calibration-in-the-large.
Calibration intercept is also a measure of calibration-in-the-large.
Moderate calibration: Calibration curves lowess based on pseudovalues
Mostrar código
# Seleccion del grupo: Stages 3-4----vdata<-dataA%>%select(id, risk2y, risk5y, time5y, eventd5y, time, eventd)%>%drop_na()primary_event<-1# A 2 años----horizon<-2# Add estimated risk and complementary log-log of it to datasetvdata$pred<-vdata$risk2ypred<-as.matrix(vdata$pred)# Calibration plot (pseudo-obs approach) ----------------------------------# First compute riskRegression::Score()score_vdata<-Score(list("csh_validation"=pred), formula =Hist(time, eventd)~1, cens.model ="km", data =vdata, conf.int =TRUE, times =horizon,# metrics = c("auc", "brier"), summary =c("ipa"), cause =primary_event, plots ="calibration")# Use pseudo-observations calculated by Score() (can alternatively use pseudo::pseudoci)pseudos<-data.frame(score_vdata$Calibration$plotframe)pseudos<-pseudos[order(pseudos$risk), ]# Use linear loess (weighted local regression with polynomial degree = 1) smoothingsmooth_pseudos<-predict(stats::loess(pseudovalue~risk, data =pseudos, degree =1, span =0.33), se =TRUE)pseudo_vals<-data.frame( obs =smooth_pseudos$fit, risk =pseudos$risk, se =smooth_pseudos$se, df =smooth_pseudos$df)%>%mutate( ll =pmax(obs-qt(0.975, df)*se, 0), ul =obs+qt(0.975, df)*se)pseudo_vals%>%ggplot(aes(x =risk, y =obs))+geom_ribbon(aes(ymin =ll, ymax =ul), fill ="grey90")+geom_abline(intercept =0, slope =1, colour ="red", linetype =2)+geom_line()+scale_y_continuous(breaks =seq(0, 1, 0.2), limits =c(0, 1))+scale_x_continuous(breaks =seq(0, 1, 0.2), limits =c(0, 1))+theme_bw()+xlab("Predicted risks")+ylab("Observed outcome proportions")->p1# A 5 años----horizon<-5# Add estimated risk and complementary log-log of it to datasetvdata$pred<-vdata$risk5ypred<-as.matrix(vdata$pred)# Calibration plot (pseudo-obs approach) ----------------------------------# First compute riskRegression::Score()score_vdata<-Score(list("csh_validation"=pred), formula =Hist(time, eventd)~1, cens.model ="km", data =vdata, conf.int =TRUE, times =horizon,# metrics = c("auc", "brier"), summary =c("ipa"), cause =primary_event, plots ="calibration")# Use pseudo-observations calculated by Score() (can alternatively use pseudo::pseudoci)pseudos<-data.frame(score_vdata$Calibration$plotframe)pseudos<-pseudos[order(pseudos$risk), ]# Use linear loess (weighted local regression with polynomial degree = 1) smoothingsmooth_pseudos<-predict(stats::loess(pseudovalue~risk, data =pseudos, degree =1, span =0.33), se =TRUE)pseudo_vals<-data.frame( obs =smooth_pseudos$fit, risk =pseudos$risk, se =smooth_pseudos$se, df =smooth_pseudos$df)%>%mutate( ll =pmax(obs-qt(0.975, df)*se, 0), ul =obs+qt(0.975, df)*se)pseudo_vals%>%ggplot(aes(x =risk, y =obs))+geom_ribbon(aes(ymin =ll, ymax =ul), fill ="grey90")+geom_abline(intercept =0, slope =1, colour ="red", linetype =2)+geom_line()+scale_y_continuous(breaks =seq(0, 1, 0.2), limits =c(0, 1))+scale_x_continuous(breaks =seq(0, 1, 0.2), limits =c(0, 1))+xlab("Predicted risks")+ylab("Observed outcome proportions")+theme_bw()->p2# Seleccion del grupo: Stages 3b-4----vdata<-dataB%>%select(id, risk2y, risk5y, time5y, eventd5y, time, eventd)%>%drop_na()primary_event<-1# A 2 años----horizon<-2# Add estimated risk and complementary log-log of it to datasetvdata$pred<-vdata$risk2ypred<-as.matrix(vdata$pred)# Calibration plot (pseudo-obs approach) ----------------------------------# First compute riskRegression::Score()score_vdata<-Score(list("csh_validation"=pred), formula =Hist(time, eventd)~1, cens.model ="km", data =vdata, conf.int =TRUE, times =horizon,# metrics = c("auc", "brier"), summary =c("ipa"), cause =primary_event, plots ="calibration")# Use pseudo-observations calculated by Score() (can alternatively use pseudo::pseudoci)pseudos<-data.frame(score_vdata$Calibration$plotframe)pseudos<-pseudos[order(pseudos$risk), ]# Use linear loess (weighted local regression with polynomial degree = 1) smoothingsmooth_pseudos<-predict(stats::loess(pseudovalue~risk, data =pseudos, degree =1, span =0.33), se =TRUE)pseudo_vals<-data.frame( obs =smooth_pseudos$fit, risk =pseudos$risk, se =smooth_pseudos$se, df =smooth_pseudos$df)%>%mutate( ll =pmax(obs-qt(0.975, df)*se, 0), ul =obs+qt(0.975, df)*se)pseudo_vals%>%ggplot(aes(x =risk, y =obs))+geom_ribbon(aes(ymin =ll, ymax =ul), fill ="grey90")+geom_abline(intercept =0, slope =1, colour ="red", linetype =2)+geom_line()+scale_y_continuous(breaks =seq(0, 1, 0.2), limits =c(0, 1))+scale_x_continuous(breaks =seq(0, 1, 0.2), limits =c(0, 1))+theme_bw()+xlab("Predicted risks")+ylab("Observed outcome proportions")->p3# A 5 años----horizon<-5# Add estimated risk and complementary log-log of it to datasetvdata$pred<-vdata$risk5ypred<-as.matrix(vdata$pred)# Calibration plot (pseudo-obs approach) ----------------------------------# First compute riskRegression::Score()score_vdata<-Score(list("csh_validation"=pred), formula =Hist(time, eventd)~1, cens.model ="km", data =vdata, conf.int =TRUE, times =horizon,# metrics = c("auc", "brier"), summary =c("ipa"), cause =primary_event, plots ="calibration")# Use pseudo-observations calculated by Score() (can alternatively use pseudo::pseudoci)pseudos<-data.frame(score_vdata$Calibration$plotframe)pseudos<-pseudos[order(pseudos$risk), ]# Use linear loess (weighted local regression with polynomial degree = 1) smoothingsmooth_pseudos<-predict(stats::loess(pseudovalue~risk, data =pseudos, degree =1, span =0.33), se =TRUE)pseudo_vals<-data.frame( obs =smooth_pseudos$fit, risk =pseudos$risk, se =smooth_pseudos$se, df =smooth_pseudos$df)%>%mutate( ll =pmax(obs-qt(0.975, df)*se, 0), ul =obs+qt(0.975, df)*se)pseudo_vals%>%ggplot(aes(x =risk, y =obs))+geom_ribbon(aes(ymin =ll, ymax =ul), fill ="grey90")+geom_abline(intercept =0, slope =1, colour ="red", linetype =2)+geom_line()+scale_y_continuous(breaks =seq(0, 1, 0.2), limits =c(0, 1))+scale_x_continuous(breaks =seq(0, 1, 0.2), limits =c(0, 1))+xlab("Predicted risks")+ylab("Observed outcome proportions")+theme_bw()->p4p1a<-p1+labs(title ="CKD Stages 3a-3b-4\n(2 year KFRE)")+theme(plot.title =element_text(hjust =0.5))p2a<-p2+labs(title ="CKD Stages 3a-3b-4\n(5 year KFRE)")+theme(plot.title =element_text(hjust =0.5))p3a<-p3+labs(title ="CKD Stages 3b-4\n(2 year KFRE)")+theme(plot.title =element_text(hjust =0.5))p4a<-p4+labs(title ="CKD Stages 3b-4\n(5 year KFRE)")+theme(plot.title =element_text(hjust =0.5))(p1a|p2a)/(p3a|p4a)+plot_annotation(tag_levels ='A')->plot_calibrationggsave(filename ="Plot_Calibration.png", device ="png", plot =plot_calibration, path =here("Figures"), scale =2, width =2100, height =2100, units ="px", dpi =300)
5.2 Fig 3
Mostrar código
plot_calibration
Figure 2: Calibration plots for each group and prediction horizon. The predicted probability is shown on the x asis, and the observed kindey failure rate is given on the y axis
5.3 Discrimination
C-index and time-dependent C/D AUC
Mostrar código
B<-1000# Seleccion del grupo: Stages 3b-4----vdata<-dataB%>%select(id, risk2y, risk5y, time5y, eventd5y, time, eventd)%>%drop_na()%>%mutate(status =factor(eventd5y, levels =c(0, 1, 2), labels =c("Alive w/o Kidney Failure", "Kidney Failure", "Death w/o Kidney Failure")), status_num =eventd5y)primary_event<-1# A 2 años----vdata$pred<-vdata$risk2yhorizon<-2pred<-as.matrix(vdata$pred)# Validation setC_vdata<-pec::cindex( object =list("CauseSpecificCox"=pred), formula =Hist(time, eventd)~1, cause =primary_event, eval.times =horizon, data =vdata)$AppCindex$CauseSpecificCox# Bootstraping C-index to calculate the bootstrap percentile confidence intervalsset.seed(777)C_boot<-function(B, data){cstat<-c()n<-nrow(data)for(iin1:B){ids<-sample(1:n, n, TRUE)data_boot<-data[ids, ]pec::cindex( object =list("CauseSpecificCox"=as.matrix(data_boot$pred)), formula =Hist(time, eventd)~1, cause =primary_event, eval.times =horizon, data =data_boot)$AppCindex$CauseSpecificCox->cstat[i]}return(cstat)}set.seed(777)C_vdata_boot<-C_boot(B =B, data =vdata)# Time-dependent AUC ---------# Validation datascore_vdata<-Score(list("csh_validation"=pred), formula =Hist(time, eventd)~1, cens.model ="km", data =vdata, conf.int =TRUE, times =horizon, metrics =c("auc"), cause =primary_event, plots ="calibration")alpha<-.05k<-3res_discr_csh<-matrix(c(## C-index# Validation CSH1C_vdata,quantile(C_vdata_boot, probs =alpha/2),quantile(C_vdata_boot, probs =1-alpha/2),## Time-dependent AUC# Validation CSH1score_vdata$AUC$score$AUC,score_vdata$AUC$score$AUC-qnorm(1-alpha/2)*score_vdata$AUC$score$se,score_vdata$AUC$score$AUC+qnorm(1-alpha/2)*score_vdata$AUC$score$se), nrow =2, ncol =3, byrow =T, dimnames =list(c("C-index, Stages 3b-4, 2 year", "Time dependent AUC, Stages 3b-4, 2 year"),c("Estimate", "Lower .95", "Upper .95")))res_discr_csh2b<-round(res_discr_csh, k)# A 5 años----vdata$pred<-vdata$risk5yhorizon<-5pred<-as.matrix(vdata$pred)# Validation setC_vdata<-pec::cindex( object =list("CauseSpecificCox"=pred), formula =Hist(time, eventd)~1, cause =primary_event, eval.times =horizon, data =vdata)$AppCindex$CauseSpecificCox# Bootstraping C-index to calculate the bootstrap percentile confidence intervalsset.seed(777)C_boot<-function(B, data){cstat<-c()n<-nrow(data)for(iin1:B){ids<-sample(1:n, n, TRUE)data_boot<-data[ids, ]pec::cindex( object =list("CauseSpecificCox"=as.matrix(data_boot$pred)), formula =Hist(time, eventd)~1, cause =primary_event, eval.times =horizon, data =data_boot)$AppCindex$CauseSpecificCox->cstat[i]}return(cstat)}set.seed(777)C_vdata_boot<-C_boot(B =B, data =vdata)# Time-dependent AUC ---------# Validation datascore_vdata<-Score(list("csh_validation"=pred), formula =Hist(time, eventd)~1, cens.model ="km", data =vdata, conf.int =TRUE, times =horizon, metrics =c("auc"), cause =primary_event, plots ="calibration")alpha<-.05k<-3res_discr_csh<-matrix(c(## C-index# Validation CSH1C_vdata,quantile(C_vdata_boot, probs =alpha/2),quantile(C_vdata_boot, probs =1-alpha/2),## Time-dependent AUC# Validation CSH1score_vdata$AUC$score$AUC,score_vdata$AUC$score$AUC-qnorm(1-alpha/2)*score_vdata$AUC$score$se,score_vdata$AUC$score$AUC+qnorm(1-alpha/2)*score_vdata$AUC$score$se), nrow =2, ncol =3, byrow =T, dimnames =list(c("C-index, Stages 3b-4, 5 year", "Time dependent AUC, Stages 3b-4, 5 year"),c("Estimate", "Lower .95", "Upper .95")))res_discr_csh5b<-round(res_discr_csh, k)
Mostrar código
# Seleccion del grupo: Stages 3a-3b-4----vdata<-dataA%>%select(id, risk2y, risk5y, time5y, eventd5y, time, eventd)%>%drop_na()%>%mutate(status =factor(eventd5y, levels =c(0, 1, 2), labels =c("Alive w/o Kidney Failure", "Kidney Failure", "Death w/o Kidney Failure")), status_num =eventd5y)primary_event<-1# A 2 años----vdata$pred<-vdata$risk2yhorizon<-2pred<-as.matrix(vdata$pred)# Validation setC_vdata<-pec::cindex( object =list("CauseSpecificCox"=pred), formula =Hist(time, eventd)~1, cause =primary_event, eval.times =horizon, data =vdata)$AppCindex$CauseSpecificCox# Bootstraping C-index to calculate the bootstrap percentile confidence intervalsset.seed(777)C_boot<-function(B, data){cstat<-c()n<-nrow(data)for(iin1:B){ids<-sample(1:n, n, TRUE)data_boot<-data[ids, ]pec::cindex( object =list("CauseSpecificCox"=as.matrix(data_boot$pred)), formula =Hist(time, eventd)~1, cause =primary_event, eval.times =horizon, data =data_boot)$AppCindex$CauseSpecificCox->cstat[i]}return(cstat)}set.seed(777)C_vdata_boot<-C_boot(B =B, data =vdata)# Time-dependent AUC ---------# Validation datascore_vdata<-Score(list("csh_validation"=pred), formula =Hist(time, eventd)~1, cens.model ="km", data =vdata, conf.int =TRUE, times =horizon, metrics =c("auc"), cause =primary_event, plots ="calibration")alpha<-.05k<-3res_discr_csh<-matrix(c(## C-index# Validation CSH1C_vdata,quantile(C_vdata_boot, probs =alpha/2),quantile(C_vdata_boot, probs =1-alpha/2),## Time-dependent AUC# Validation CSH1score_vdata$AUC$score$AUC,score_vdata$AUC$score$AUC-qnorm(1-alpha/2)*score_vdata$AUC$score$se,score_vdata$AUC$score$AUC+qnorm(1-alpha/2)*score_vdata$AUC$score$se), nrow =2, ncol =3, byrow =T, dimnames =list(c("C-index, Stages 3-4, 2 year", "Time dependent AUC, Stages 3-4, 2 year"),c("Estimate", "Lower .95", "Upper .95")))res_discr_csh2a<-round(res_discr_csh, k)# A 5 años----vdata$pred<-vdata$risk5yhorizon<-5pred<-as.matrix(vdata$pred)# Validation setC_vdata<-pec::cindex( object =list("CauseSpecificCox"=pred), formula =Hist(time, eventd)~1, cause =primary_event, eval.times =horizon, data =vdata)$AppCindex$CauseSpecificCox# Bootstraping C-index to calculate the bootstrap percentile confidence intervalsset.seed(777)C_boot<-function(B, data){cstat<-c()n<-nrow(data)for(iin1:B){ids<-sample(1:n, n, TRUE)data_boot<-data[ids, ]pec::cindex( object =list("CauseSpecificCox"=as.matrix(data_boot$pred)), formula =Hist(time, eventd)~1, cause =primary_event, eval.times =horizon, data =data_boot)$AppCindex$CauseSpecificCox->cstat[i]}return(cstat)}set.seed(777)C_vdata_boot<-C_boot(B =B, data =vdata)# Time-dependent AUC ---------# Validation datascore_vdata<-Score(list("csh_validation"=pred), formula =Hist(time, eventd)~1, cens.model ="km", data =vdata, conf.int =TRUE, times =horizon, metrics =c("auc"), cause =primary_event, plots ="calibration")alpha<-.05k<-3res_discr_csh<-matrix(c(## C-index# Validation CSH1C_vdata,quantile(C_vdata_boot, probs =alpha/2),quantile(C_vdata_boot, probs =1-alpha/2),## Time-dependent AUC# Validation CSH1score_vdata$AUC$score$AUC,score_vdata$AUC$score$AUC-qnorm(1-alpha/2)*score_vdata$AUC$score$se,score_vdata$AUC$score$AUC+qnorm(1-alpha/2)*score_vdata$AUC$score$se), nrow =2, ncol =3, byrow =T, dimnames =list(c("C-index, Stages 3-4, 5 year", "Time dependent AUC, Stages 3-4, 5 year"),c("Estimate", "Lower .95", "Upper .95")))res_discr_csh5a<-round(res_discr_csh, k)
Mostrar código
# Average predicted riskavg_pred<-cbind("metrica"="Average predicted risk", avg_pred2a, avg_pred5a, avg_pred2b, avg_pred5b)colnames(avg_pred)<-c("metrica", "est2ya", "est5ya", "est2yb", "est5yb")avg_pred%>%as_tibble()%>%mutate( est2ya =as.character(glue("{est2ya}%")), est5ya =as.character(glue("{est5ya}%")), est2yb =as.character(glue("{est2yb}%")), est5yb =as.character(glue("{est5yb}%")))->avg_pred; avg_pred# Average observed proportionavg_obs<-cbind("metrica"="Average observed proportion (95% CI)", avg_obs2a, avg_obs5a, avg_obs2b, avg_obs5b)colnames(avg_obs)<-c("metrica", "OE2a", "ll2a", "ul2a", "OE5a", "ll5a", "ul5a", "OE2b", "ll2b", "ul2b", "OE5b", "ll5b", "ul5b")avg_obs%>%as_tibble()%>%mutate( est2ya =as.character(glue("{OE2a}% ({ll2a}% to {ul2a}%)")), est5ya =as.character(glue("{OE5a}% ({ll5a}% to {ul5a}%)")), est2yb =as.character(glue("{OE2b}% ({ll2b}% to {ul2b}%)")), est5yb =as.character(glue("{OE5b}% ({ll5b}% to {ul5b}%)")))%>%select(metrica, starts_with("est"))->avg_obs; avg_obs# OE summaryOE_summary<-cbind("metrica"="O/E ratio (95% CI)", OE_summary2a, OE_summary5a, OE_summary2b, OE_summary5b); OE_summarycolnames(OE_summary)<-c("metrica", "OE2a", "ll2a", "ul2a", "OE5a", "ll5a", "ul5a", "OE2b", "ll2b", "ul2b", "OE5b", "ll5b", "ul5b")OE_summary%>%as_tibble()%>%mutate( est2ya =as.character(glue("{OE2a} ({ll2a} to {ul2a})")), est5ya =as.character(glue("{OE5a} ({ll5a} to {ul5a})")), est2yb =as.character(glue("{OE2b} ({ll2b} to {ul2b})")), est5yb =as.character(glue("{OE5b} ({ll5b} to {ul5b})")))%>%select(metrica, starts_with("est"))->OE_summary; OE_summary# Calibration slope e interceptres_cal<-cbind("metrica"=c("Calibration intercept (95% CI)", " Calibration slope (95% CI)"), res_cal2a, res_cal5a, res_cal2b, res_cal5b); res_calcolnames(res_cal)<-c("metrica", "OE2a", "ll2a", "ul2a", "OE5a", "ll5a", "ul5a", "OE2b", "ll2b", "ul2b", "OE5b", "ll5b", "ul5b")res_cal<-res_cal%>%as_tibble()%>%mutate( est2ya =as.character(glue("{OE2a} ({ll2a} to {ul2a})")), est5ya =as.character(glue("{OE5a} ({ll5a} to {ul5a})")), est2yb =as.character(glue("{OE2b} ({ll2b} to {ul2b})")), est5yb =as.character(glue("{OE5b} ({ll5b} to {ul5b})")))%>%select(metrica, starts_with("est"))rbind(avg_pred, avg_obs, OE_summary, res_cal)->table_performance; table_performance# C-index y C/D AUCres_discr_csh<-cbind("metrica"=c("C-index up to t-years (95% CI)", " C/D AUC, at t years (95% CI)"), res_discr_csh2a, res_discr_csh5a, res_discr_csh2b, res_discr_csh5b); res_discr_cshcolnames(res_discr_csh)<-c("metrica", "OE2a", "ll2a", "ul2a", "OE5a", "ll5a", "ul5a", "OE2b", "ll2b", "ul2b", "OE5b", "ll5b", "ul5b")res_discr_csh<-res_discr_csh%>%as_tibble()%>%mutate( est2ya =as.character(glue("{OE2a} ({ll2a} to {ul2a})")), est5ya =as.character(glue("{OE5a} ({ll5a} to {ul5a})")), est2yb =as.character(glue("{OE2b} ({ll2b} to {ul2b})")), est5yb =as.character(glue("{OE5b} ({ll5b} to {ul5b})")))%>%select(metrica, starts_with("est"))rbind(avg_pred, avg_obs, OE_summary, res_cal, res_discr_csh)->table_performance; table_performance
Mostrar código
table_performance%>%mutate(grupo =c(rep("Calibration", 5), rep("Discrimination", 2)))%>%relocate(grupo, .before ="metrica")%>%bind_rows()%>%as_grouped_data(groups ="grupo")%>%flextable::as_flextable(hide_grouplabel =TRUE)%>%set_header_labels( metrica ="Validation aspect and performance measure", est2ya ="t = 2 year", est5ya ="t = 5 year", est2yb ="t = 2 year", est5yb ="t = 5 year")%>%add_header_row( values =c("Validation aspect and performance measure", "CKD Stages 3a-3b-4", "CKD Stages 3b-4"), colwidths =c(1, 2, 2))%>%merge_v(j =1, part ="header")%>%bold(i =c(1, 7))%>%autofit()%>%set_caption("Table 2. Performance measures of KFRE in the external dataset of patients with CKD Stages 3a-4 and 3b-4")%>%theme_booktabs()%>%bold(bold =TRUE, part ="header")->table_perf_finaltable_perf_final%>%flextable::save_as_docx(path =here("Tables/Table2.docx"))
5.4 Table 2
Mostrar código
table_perf_final
Validation aspect and performance measure
CKD Stages 3a-3b-4
CKD Stages 3b-4
t = 2 year
t = 5 year
t = 2 year
t = 5 year
Calibration
Average predicted risk
0.96%
3.18%
2.36%
7.66%
Average observed proportion (95% CI)
1.52% (1.24% to 1.79%)
3.37% (2.95% to 3.8%)
3.15% (2.5% to 3.79%)
6.86% (5.89% to 7.84%)
O/E ratio (95% CI)
1.57 (1.39 to 1.76)
1.06 (0.93 to 1.19)
1.33 (1.13 to 1.54)
0.9 (0.75 to 1.04)
Calibration intercept (95% CI)
0.18 (-0.1 to 0.45)
-0.26 (-0.45 to -0.07)
0.16 (-0.12 to 0.44)
-0.29 (-0.48 to -0.1)
Calibration slope (95% CI)
0.79 (0.61 to 0.96)
0.75 (0.65 to 0.86)
0.82 (0.6 to 1.03)
0.79 (0.66 to 0.92)
Discrimination
C-index up to t-years (95% CI)
0.853 (0.812 to 0.892)
0.845 (0.818 to 0.872)
0.848 (0.803 to 0.885)
0.827 (0.796 to 0.857)
C/D AUC, at t years (95% CI)
0.855 (0.816 to 0.895)
0.847 (0.819 to 0.875)
0.853 (0.811 to 0.895)
0.836 (0.803 to 0.869)
6 Supplementary tables
6.1 Table S1
Mostrar código
table_kfre<-data.frame( pred =c("2-years", "5-years"), eq =c("$1-{0.9832}^{e^{(-0.2201\times(\frac{age}{10}-7.036)+0.2467\times(male-0.5642)-0.5567\times(\frac{eGFR}{5}-7.222)+0.4510\times(log{(ACR)}-5.137))}}$", "$1-{0.9365}^{e^{(-0.2201\times(\frac{age}{10}-7.036)+0.2467\times(male-0.5642)-0.5567\times(\frac{eGFR}{5}-7.222)+0.4510\times(log{(ACR)}-5.137))}}$"))
Mostrar código
knitr::kable(table_kfre, escape =TRUE, col.names =c("Prediction horizons", "Original regional equation calibrated for predicted risk of kidney failure"), caption ="**Table S1.** KFRE equations externally validated by the study")
**Table S1.** KFRE equations externally validated by the study
Prediction horizons
Original regional equation calibrated for predicted risk of kidney failure
table_coding<-data.frame( Variable =c("age", "male", "eGFR_ckdepi", "acr"), Coding =c("integer number that indicates the age in completed years", "1 = male; 0 = female", "estimated glomerular filtration rate obtained by CKD-EPI formula in $ml/min/1.73m^2$", "albumin-to-creatinine ratio in mg/g"))
Mostrar código
knitr::kable(table_coding, escape =TRUE, caption ="*Table S2.* Coding of variables")
*Table S2.* Coding of variables
Variable
Coding
age
integer number that indicates the age in completed years
male
1 = male; 0 = female
eGFR_ckdepi
estimated glomerular filtration rate obtained by CKD-EPI formula in $ml/min/1.73m^2$
Table S6. Distribution of CKD 3a-4 patients included in the analysis within 17 health facilities of the EsSalud Rebagliati Network
CKD Stages 3a-3b-4
Overall
Outcome
Characteristic
N = 7,519
Alive w/o Kidney Failure, N = 5,880
Death w/o Kidney Failure, N = 1,400
Kidney Failure, N = 239
Healthcare center
C. M. MALA
74 (100.0%)
48 (64.9%)
18 (24.3%)
8 (10.8%)
CAP II LURIN
58 (100.0%)
52 (89.7%)
6 (10.3%)
0 (0.0%)
CAP III LOS PROCEDES DE SJM
403 (100.0%)
319 (79.2%)
73 (18.1%)
11 (2.7%)
CAP III SAN ISIDRO
287 (100.0%)
204 (71.1%)
71 (24.7%)
12 (4.2%)
CAP III SURQUILLO
308 (100.0%)
242 (78.6%)
53 (17.2%)
13 (4.2%)
HOSP. II ANGAMOS
1,471 (100.0%)
1,088 (74.0%)
298 (20.3%)
85 (5.8%)
HOSPITAL EDGARDO REBAGLIATI MARTINS
34 (100.0%)
20 (58.8%)
11 (32.4%)
3 (8.8%)
HOSPITAL I CARLOS ALCANTARA B.
1,032 (100.0%)
858 (83.1%)
160 (15.5%)
14 (1.4%)
HOSPITAL I ULDARICO ROCCA F.
207 (100.0%)
158 (76.3%)
42 (20.3%)
7 (3.4%)
HOSPITAL II CANETE
258 (100.0%)
179 (69.4%)
73 (28.3%)
6 (2.3%)
POLICLINICO CHEQUEOS LARCO
77 (100.0%)
71 (92.2%)
5 (6.5%)
1 (1.3%)
POLICLINICO CHINCHA
1,289 (100.0%)
1,058 (82.1%)
202 (15.7%)
29 (2.2%)
POLICLINICO JUAN JOSE RODRIGUEZ LAZO
679 (100.0%)
497 (73.2%)
166 (24.4%)
16 (2.4%)
POLICLINICO PABLO BERMUDEZ
752 (100.0%)
586 (77.9%)
143 (19.0%)
23 (3.1%)
POLICLINICO PROCERES
454 (100.0%)
381 (83.9%)
64 (14.1%)
9 (2.0%)
POLICLINICO SANTA CRUZ
102 (100.0%)
90 (88.2%)
10 (9.8%)
2 (2.0%)
POLICLINICO VILLA MARIA
34 (100.0%)
29 (85.3%)
5 (14.7%)
0 (0.0%)
1 n (%)
6.7 Table S7
Mostrar código
tbl_merge(list(tab1b, tab2b), tab_spanner =c("Overall", "Outcome"))%>%modify_caption("Table S7. Distribution of CKD 3b-4 patients included in the analysis within 17 health facilities of the EsSalud Rebagliati Network")->tabBtabB%>%as_flex_table()%>%theme_booktabs()%>%bold(bold =TRUE, part ="header")%>%flextable::save_as_docx(path =here("Tables/TableS7.docx"))
Table S7. Distribution of CKD 3b-4 patients included in the analysis within 17 health facilities of the EsSalud Rebagliati Network
CKD Stages 3b-4
Overall
Outcome
Characteristic
N = 2,798
Alive w/o Kidney Failure, N = 1,933
Death w/o Kidney Failure, N = 683
Kidney Failure, N = 182
Healthcare center
C. M. MALA
26 (100.0%)
10 (38.5%)
10 (38.5%)
6 (23.1%)
CAP II LURIN
12 (100.0%)
8 (66.7%)
4 (33.3%)
0 (0.0%)
CAP III LOS PROCEDES DE SJM
129 (100.0%)
87 (67.4%)
32 (24.8%)
10 (7.8%)
CAP III SAN ISIDRO
113 (100.0%)
68 (60.2%)
35 (31.0%)
10 (8.8%)
CAP III SURQUILLO
115 (100.0%)
81 (70.4%)
24 (20.9%)
10 (8.7%)
HOSP. II ANGAMOS
921 (100.0%)
644 (69.9%)
212 (23.0%)
65 (7.1%)
HOSPITAL EDGARDO REBAGLIATI MARTINS
22 (100.0%)
11 (50.0%)
8 (36.4%)
3 (13.6%)
HOSPITAL I CARLOS ALCANTARA B.
271 (100.0%)
194 (71.6%)
67 (24.7%)
10 (3.7%)
HOSPITAL I ULDARICO ROCCA F.
70 (100.0%)
47 (67.1%)
17 (24.3%)
6 (8.6%)
HOSPITAL II CANETE
90 (100.0%)
59 (65.6%)
28 (31.1%)
3 (3.3%)
POLICLINICO CHEQUEOS LARCO
15 (100.0%)
13 (86.7%)
1 (6.7%)
1 (6.7%)
POLICLINICO CHINCHA
408 (100.0%)
294 (72.1%)
93 (22.8%)
21 (5.1%)
POLICLINICO JUAN JOSE RODRIGUEZ LAZO
216 (100.0%)
134 (62.0%)
70 (32.4%)
12 (5.6%)
POLICLINICO PABLO BERMUDEZ
215 (100.0%)
152 (70.7%)
47 (21.9%)
16 (7.4%)
POLICLINICO PROCERES
144 (100.0%)
106 (73.6%)
30 (20.8%)
8 (5.6%)
POLICLINICO SANTA CRUZ
23 (100.0%)
18 (78.3%)
4 (17.4%)
1 (4.3%)
POLICLINICO VILLA MARIA
8 (100.0%)
7 (87.5%)
1 (12.5%)
0 (0.0%)
1 n (%)
6.8 Table S8
Mostrar código
# Selection of group of patients----group<-c("Stages 3-4")# Creacion de datsets para IA a 5 y 2 añosdata%>%filter(ckd_stage==group)->data_filtvdata.w<-crprep( Tstop ="time", status ="eventd", trans =c(1, 2), id ="id", keep =c("age", "male", "eGFR_mdrd", "acr", "risk2y", "risk5y"), data =data_filt)vdata.w1<-vdata.w%>%filter(failcode==1)vdata.w2<-vdata.w%>%filter(failcode==2)# For kidney failuremfit_vdata1<-survfit(Surv(Tstart, Tstop, status==1)~1, data =vdata.w1, weights =weight.cens)smfit_vdata1<-summary(mfit_vdata1, times =c(1, 2, 3, 4, 5))res_ci_stg1<-cbind(100*(1-smfit_vdata1$surv),100*(1-smfit_vdata1$upper),100*(1-smfit_vdata1$lower))res_ci_stg1<-round(res_ci_stg1, 2)rownames(res_ci_stg1)<-c("1-year", "2-year","3-year", "4-year","5-year")colnames(res_ci_stg1)<-c("Estimate", "Lower .95","Upper .95")# For death without kidney failure mfit_vdata2<-survfit(Surv(Tstart, Tstop, status==2)~1, data =vdata.w2, weights =weight.cens)smfit_vdata2<-summary(mfit_vdata2, times =c(1, 2, 3, 4, 5))res_ci_stg2<-cbind(100*(1-smfit_vdata2$surv),100*(1-smfit_vdata2$upper),100*(1-smfit_vdata2$lower))res_ci_stg2<-round(res_ci_stg2, 2)rownames(res_ci_stg2)<-c("1-year", "2-year","3-year", "4-year","5-year")colnames(res_ci_stg2)<-c("Estimate", "Lower .95","Upper .95")# Selection of group of patients----group<-c("Stages 3b-4")# Creation of datasets for Cumulative Incidence at 2- and 5-yearsadata%>%filter(ckd_stage2==group)->data_filtvdata.w<-crprep( Tstop ="time5y", status ="eventd5y", trans =c(1, 2), id ="id", keep =c("age", "male", "eGFR_mdrd", "acr", "risk2y", "risk5y"), data =data_filt)vdata.w1<-vdata.w%>%filter(failcode==1)vdata.w2<-vdata.w%>%filter(failcode==2)# For kidney failuremfit_vdata3<-survfit(Surv(Tstart, Tstop, status==1)~1, data =vdata.w1, weights =weight.cens)smfit_vdata3<-summary(mfit_vdata3, times =c(1, 2, 3, 4, 5))res_ci_stg3<-cbind(100*(1-smfit_vdata3$surv),100*(1-smfit_vdata3$upper),100*(1-smfit_vdata3$lower))res_ci_stg3<-round(res_ci_stg3, 2)rownames(res_ci_stg3)<-c("1-year", "2-year","3-year", "4-year","5-year")colnames(res_ci_stg3)<-c("Estimate", "Lower .95","Upper .95")# For death without kidney failure mfit_vdata4<-survfit(Surv(Tstart, Tstop, status==2)~1, data =vdata.w2, weights =weight.cens)smfit_vdata4<-summary(mfit_vdata4, times =c(1, 2, 3, 4, 5))res_ci_stg4<-cbind(100*(1-smfit_vdata4$surv),100*(1-smfit_vdata4$upper),100*(1-smfit_vdata4$lower))res_ci_stg4<-round(res_ci_stg4, 2)rownames(res_ci_stg4)<-c("1-year", "2-year","3-year", "4-year","5-year")colnames(res_ci_stg4)<-c("Estimate", "Lower .95","Upper .95")# Table for 3a-3b-4 CKD patients----res_ci<-cbind(res_ci_stg1, res_ci_stg2)res_ci%>%as_tibble(rownames ="Year")%>%mutate( est1 =str_glue("{Estimate}%"), est1.ci =str_glue(" ({`Lower .95`}% to {`Upper .95`}%)"), est2 =str_glue("{V4}%"), est2.ci =str_glue(" ({V5}% to {V6}%)"))%>%select(Year, est1, est1.ci, est2, est2.ci)%>%mutate(across(everything(), ~as.character(.x)))->res_ci2Ares_ci2A%>%flextable()%>%add_header_row(values =c("Year", "Kidney failure", "Death without kidney failure"), colwidths =c(1, 2, 2))%>%set_header_labels( est1 ="%", est1.ci ="95% CI", est2 ="%", est2.ci ="95% CI")%>%merge_v(j =1, part ="header")%>%set_caption("Table S8. Cumulative incidence of kidney failure and death without kidney failure in patients with CKD stages 3a-3b-4")%>%theme_box()%>%autofit()->tab_cif1tab_cif1%>%save_as_docx(path =here("Tables/TableS8.docx"))
Mostrar código
kable(res_ci2A, caption ="Table S8. Cumulative incidence of Kidney Failure and Death w/o Kidney Failure in patients with CKD Stages 3a-3b-4", col.names =c("", "%", "95%CI", "%", "95%CI"))%>%kable_styling("striped", position ="center")%>%add_header_above(c(" "=1, "Kidney Failure"=2, "Death w/o Kidney Failure"=2))
Table S8. Cumulative incidence of Kidney Failure and Death w/o Kidney Failure in patients with CKD Stages 3a-3b-4
Kidney Failure
Death w/o Kidney Failure
%
95%CI
%
95%CI
1-year
0.68%
(0.49% to 0.86%)
3.82%
(3.38% to 4.25%)
2-year
1.52%
(1.24% to 1.79%)
7.49%
(6.89% to 8.08%)
3-year
2.23%
(1.9% to 2.57%)
11.2%
(10.48% to 11.91%)
4-year
2.88%
(2.5% to 3.26%)
15.58%
(14.75% to 16.41%)
5-year
3.37%
(2.95% to 3.8%)
20.46%
(19.48% to 21.42%)
The estimated cumulative incidence of renal failure at 2 and 5 years was 1.52% ( (1.24% to 1.79%)) and 3.37% ( (2.95% to 3.8%)) in patients with stage 3a, 3b or 4 CKD, respectively (Table S8).
6.9 Tabla S9
Mostrar código
# Tabla for 3b-4----res_ci<-cbind(res_ci_stg3, res_ci_stg4)res_ci%>%as_tibble(rownames ="Year")%>%mutate( est1 =str_glue("{Estimate}%"), est1.ci =str_glue(" ({`Lower .95`}% to {`Upper .95`}%)"), est2 =str_glue("{V4}%"), est2.ci =str_glue(" ({V5}% to {V6}%)"))%>%select(Year, est1, est1.ci, est2, est2.ci)%>%mutate(across(everything(), ~as.character(.x)))->res_ci2Bres_ci2B%>%flextable()%>%add_header_row(values =c("Year", "Kidney failure", "Death without kidney failure"), colwidths =c(1, 2, 2))%>%set_header_labels( est1 ="%", est1.ci ="95% CI", est2 ="%", est2.ci ="95% CI")%>%merge_v(j =1, part ="header")%>%set_caption("Table S9. Cumulative incidence of kidney failure and death without kidney failure in patients with CKD stages 3b-4")%>%theme_box()%>%autofit()->tab_cif2tab_cif2%>%save_as_docx(path =here("Tables/TableS9.docx"))
Mostrar código
kable(res_ci2B, caption ="Table S9. Cumulative incidence of Kidney Failure and Death w/o Kidney Failure in patients with CKD Stages 3b-4", col.names =c("", "%", "95%CI", "%", "95%CI"))%>%kable_styling("striped", position ="center")%>%add_header_above(c(" "=1, "Kidney Failure"=2, "Death w/o Kidney Failure"=2))
Table S9. Cumulative incidence of Kidney Failure and Death w/o Kidney Failure in patients with CKD Stages 3b-4
Kidney Failure
Death w/o Kidney Failure
%
95%CI
%
95%CI
1-year
1.54%
(1.08% to 1.99%)
6%
(5.12% to 6.88%)
2-year
3.15%
(2.5% to 3.79%)
10.72%
(9.57% to 11.86%)
3-year
4.72%
(3.93% to 5.5%)
15.51%
(14.16% to 16.84%)
4-year
6.02%
(5.12% to 6.9%)
21.03%
(19.48% to 22.55%)
5-year
6.86%
(5.89% to 7.83%)
26.59%
(24.84% to 28.3%)
In patients with stage 3b-4 CKD, the cumulative incidences at 2 and 5 years were 10.72% ( (9.57% to 11.86%)) and 26.59% ( (24.84% to 28.3%)), respectively. Table S9 details the cumulative incidences of renal failure from 1 to 5 years in both cohorts.
6.10 Fig S1
Mostrar código
data%>%filter(ckd_stage=="Stages 3-4")%>%mutate(evento ="CKD 3-4")->dataAdata%>%filter(ckd_stage2=="Stages 3b-4")%>%mutate(evento ="CKD 3b-4")->dataBdataA%>%bind_rows(dataB)->data2data2%>%ggplot(aes(x =evento, y =age, fill =evento))+sm_raincloud()+labs(x ="", y ="Age, years")+coord_cartesian(ylim =c(0, 100))->p1data2%>%ggplot(aes(x =evento, y =eGFR_ckdepi, fill =evento))+sm_raincloud()+labs(x ="", y ="eGFR, ml/min/1.73m<sup>2</sup>")+coord_cartesian(ylim =c(0, 60))+theme(axis.title.y =element_markdown())->p2data2%>%ggplot(aes(x =evento, y =acr, fill =evento))+sm_raincloud()+labs(x ="", y ="ACR, mg/g")+coord_cartesian(ylim =c(0, 150000))->p3data2%>%ggplot(aes(x =evento, y =acr, fill =evento))+sm_raincloud()+scale_y_continuous(trans ="log10", labels =label_log())+labs(x ="", y ="ACR, mg/g (log scale)")->p4(p1+p2)/(p3+p4)+plot_annotation(tag_levels ="A")->plot_dist_varsggsave(filename ="plot_dist_vars.png", plot =plot_dist_vars, device ="png", path =here("Figures/"), scale =2, width =18, height =18, units ="cm", dpi =900)
Figure 3: Distribution of four variables of KFRE equation in CKD 3-4 and CKD 3b-4 patientes. (A) age in years, (B) estimate glomerural filtration rate (eGFR) according to CKD-EPI formula, (C) urine albumin to creatinine ratio (ACR) expressed in original scale and in (D) natural logarithm scale for better comparison of the distributions.
6.11 Fig S2
Mostrar código
data%>%filter(ckd_stage=="Stages 3-4")%>%mutate( evento2y =case_when(eventd2ylab%in%c("Alive w/o Kidney Failure", "Death w/o Kidney Failure")~"No", eventd2ylab=="Kidney Failure"~"Yes", TRUE~as.character(NA)), evento5y =case_when(eventd5ylab%in%c("Alive w/o Kidney Failure", "Death w/o Kidney Failure")~"No", eventd5ylab=="Kidney Failure"~"Yes", TRUE~as.character(NA)))->dataAdata%>%filter(ckd_stage2=="Stages 3b-4")%>%mutate( evento2y =case_when(eventd2ylab%in%c("Alive w/o Kidney Failure", "Death w/o Kidney Failure")~"No", eventd2ylab=="Kidney Failure"~"Yes", TRUE~as.character(NA)), evento5y =case_when(eventd5ylab%in%c("Alive w/o Kidney Failure", "Death w/o Kidney Failure")~"No", eventd5ylab=="Kidney Failure"~"Yes", TRUE~as.character(NA)))->dataBdataA%>%ggplot(aes(x =risk2y, color =evento2y, fill =evento2y))+geom_histogram(alpha =0.4)+labs(x ="2-years predicted risk by KFRE", y ="Frequency count", color ="Kidney failure", fill ="Kidney failure", title ="CKD stages 3-4")+theme_classic()+theme(plot.title =element_text(hjust =0.5, face ="bold"))+coord_cartesian(xlim =c(0, 1))->p1dataA%>%ggplot(aes(x =risk5y, color =evento5y, fill =evento5y))+geom_histogram(alpha =0.4)+labs(x ="5-years predicted risk by KFRE", y ="Frequency count", color ="Kidney failure", fill ="Kidney failure")+theme_classic()+coord_cartesian(xlim =c(0, 1))->p2dataB%>%ggplot(aes(x =risk2y, color =evento2y, fill =evento2y))+geom_histogram(alpha =0.4)+labs(x ="2-years predicted risk by KFRE", y ="Frequency count", color ="Kidney failure", fill ="Kidney failure", title ="CKD stages 3b-4")+theme_classic()+theme(plot.title =element_text(hjust =0.5, face ="bold"))+coord_cartesian(xlim =c(0, 1))->p3dataB%>%ggplot(aes(x =risk5y, color =evento5y, fill =evento5y))+geom_histogram(alpha =0.4)+labs(x ="5-years predicted risk by KFRE", y ="Frequency count", color ="Kidney failure", fill ="Kidney failure")+theme_classic()+coord_cartesian(xlim =c(0, 1))->p4((p1/p2)|(p3/p4))+plot_layout(guides ="collect")+plot_annotation(tag_levels ="A")->plot_dist_risksggsave(filename ="plot_dist_risks.png", plot =plot_dist_risks, device ="png", path =here("Figures/"), scale =2, width =18, height =18, units ="cm", dpi =900)
